---
title: "Table 2 and tables in Online appendix IV"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Note

The code in the files "ecolRxC_accuracy_IPF.html", "ecolRxC_accuracy_Thomsen_logit_no_aprox.html", "ecolRxC_accuracy_Thomsen_logit_Yule_aprox.html",  "ecolRxC_accuracy_Thomsen_probit_no_aprox.html" and "ecolRxC_accuracy_Thomsen_probit_Yule_aprox.html" must be previously run for this code to work. The package **tidyverse** is needed for this code to run.

## Introduction

This document contains the code to replicate Table 2 in the paper and the figures in the tables in Online appendix IV.

### Reproducible code

#### Data
```{r}
library(tidyverse)
datos1 <- read.csv("summary_ecolRxC_solutions_IPF.csv", sep =";", stringsAsFactors = F)
datos2 <- read.csv("summary_ecolRxC_solutions_Thomsen_logit_FALSE.csv", sep =";", 
                   stringsAsFactors = F)
datos3 <- read.csv("summary_ecolRxC_solutions_Thomsen_logit_TRUE.csv", sep =";", stringsAsFactors = F)
datos4 <- read.csv("summary_ecolRxC_solutions_Thomsen_probit_FALSE.csv", sep =";",
                  stringsAsFactors = F)
datos5 <- read.csv("summary_ecolRxC_solutions_Thomsen_probit_TRUE.csv", sep =";", stringsAsFactors = F)
```

#### EI.logit (Top panel Table 2)
```{r}
datos1.1 <- datos1 %>% select(year, EI.ecolRxC.logit, EI.ecolRxC.logit.global,
                             EI.ecolRxC.logit.Yule, EI.ecolRxC.logit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos2.1 <- datos2 %>% select(EI.ecolRxC.AVCR, EI.ecolRxC.Average) 
colnames(datos2.1) <-c("ecolRxC", "ecol-biN")

datos3.1 <- datos3 %>% select(EI.ecolRxC.AVCR, EI.ecolRxC.Average) 
colnames(datos3.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos2.1, datos3.1)

EI.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

orden <- order(as.matrix(NZ.SCO[1, ]), decreasing = TRUE)
orden <- c(1, orden + 1)

NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EI.errors <- rbind(EI.errors, NZ.SCO)[, orden]
write.table(EI.errors, "EI_logit.csv", sep = ";", row.names = FALSE)
```

#### EI.probit (Bottom panel Table 2)
```{r}
datos1.1 <- datos1 %>% select(year, EI.ecolRxC.probit, EI.ecolRxC.probit.global,
                             EI.ecolRxC.probit.Yule, EI.ecolRxC.probit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos4.1 <- datos4 %>% select(EI.ecolRxC.AVCR, EI.ecolRxC.Average) 
colnames(datos4.1) <-c("ecolRxC", "ecol-biN")

datos5.1 <- datos5 %>% select(EI.ecolRxC.AVCR, EI.ecolRxC.Average) 
colnames(datos5.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos4.1, datos5.1)

EI.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EI.errors <- rbind(EI.errors, NZ.SCO)[, orden]

write.table(EI.errors, "EI_probit.csv", sep = ";", row.names = FALSE)
```



#### EPW.logit (Appendix IV)
```{r}
datos1.1 <- datos1 %>% select(year, EPW.ecolRxC.logit, EPW.ecolRxC.logit.global,
                             EPW.ecolRxC.logit.Yule, EPW.ecolRxC.logit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos2.1 <- datos2 %>% select(EPW.ecolRxC.AVCR, EPW.ecolRxC.Average) 
colnames(datos2.1) <-c("ecolRxC", "ecol-biN")

datos3.1 <- datos3 %>% select(EPW.ecolRxC.AVCR, EPW.ecolRxC.Average) 
colnames(datos3.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos2.1, datos3.1)

EPW.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EPW.errors <- rbind(EPW.errors, NZ.SCO)[, orden]
write.table(EPW.errors, "EPW_logit.csv", sep = ";", row.names = FALSE)
```

#### EPW.probit (Appendix IV)
```{r}
datos1.1 <- datos1 %>% select(year, EPW.ecolRxC.probit, EPW.ecolRxC.probit.global,
                             EPW.ecolRxC.probit.Yule, EPW.ecolRxC.probit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos4.1 <- datos4 %>% select(EPW.ecolRxC.AVCR, EPW.ecolRxC.Average) 
colnames(datos4.1) <-c("ecolRxC", "ecol-biN")

datos5.1 <- datos5 %>% select(EPW.ecolRxC.AVCR, EPW.ecolRxC.Average) 
colnames(datos5.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos4.1, datos5.1)


EPW.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EPW.errors <- rbind(EPW.errors, NZ.SCO)[, orden]

write.table(EPW.errors, "EPW_probit.csv", sep = ";", row.names = FALSE)
```



#### EQ.logit (Appendix IV)
```{r}
datos1.1 <- datos1 %>% select(year, EQ.ecolRxC.logit, EQ.ecolRxC.logit.global,
                             EQ.ecolRxC.logit.Yule, EQ.ecolRxC.logit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos2.1 <- datos2 %>% select(EQ.ecolRxC.AVCR, EQ.ecolRxC.Average) 
colnames(datos2.1) <-c("ecolRxC", "ecol-biN")

datos3.1 <- datos3 %>% select(EQ.ecolRxC.AVCR, EQ.ecolRxC.Average) 
colnames(datos3.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos2.1, datos3.1)

EQ.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EQ.errors <- rbind(EQ.errors, NZ.SCO)[, orden]
write.table(EQ.errors, "EQ_logit.csv", sep = ";", row.names = FALSE)
```

#### EQ.probit (Appendix IV)
```{r}
datos1.1 <- datos1 %>% select(year, EQ.ecolRxC.probit, EQ.ecolRxC.probit.global,
                             EQ.ecolRxC.probit.Yule, EQ.ecolRxC.probit.global.Yule) 
colnames(datos1.1) <-c("year", "VTR-local", "VTR", "VTR-local-Yule", "VTR-Yule")

datos4.1 <- datos4 %>% select(EQ.ecolRxC.AVCR, EQ.ecolRxC.Average) 
colnames(datos4.1) <-c("ecolRxC", "ecol-biN")

datos5.1 <- datos5 %>% select(EQ.ecolRxC.AVCR, EQ.ecolRxC.Average) 
colnames(datos5.1) <-c("ecolRxC-Yule", "ecol")

datos <- cbind(datos1.1, datos4.1, datos5.1)

EQ.errors <- datos %>% group_by(year) %>% 
  summarise(across(everything(), mean, na.rm = TRUE))

NZ.SCO <-  datos %>% select(`VTR-local`:ecol) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
NZ.SCO <- cbind(year = "NZ+SCO", NZ.SCO)

EQ.errors <- rbind(EQ.errors, NZ.SCO)[, orden]
write.table(EQ.errors, "EQ_probit.csv", sep = ";", row.names = FALSE)
```
